Influence of pions and hyperons on stellar black hole formation 
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We present numerical simulations of stellar core-collapse with spherically symmetric, general 
relativistic hydrodynamics up to black hole formation. Using the CoCoNuT code, with a newly 
developed grey leakage scheme for the neutrino treatment, we investigate the effects of including 
pions and A-hyperons into the equation of state at high densities and temperatures on the black 
hole formation process. Results show non-negligible differences between the models with reference 
equation of state without any additional particles and models with the extended ones. For the 
latter, the maximum masses supported by the proto-neutron star are smaller and the collapse to a 
black hole occurs earlier. A phase transition to hyperonic matter is observed when the progenitor 
allows for a high enough accretion rate onto the proto-neutron star. Rough estimates of neutrino 
luminosity from these collapses are given, too. 

PACS numbers: 97.60.Lf, 26.50,+x, 97.60.Bw 



I. INTRODUCTION 

Supernovae and hypernovae are among the most spec- 
tacular events in the observable universe with an enor- 
mous amount of energy involved. Numerous studies have 
been undertaken to understand their mechanisms, for re- 
cent reviews, see e.g. [H-Q and references therein. At the 
origin of an iron core-collapse supernovae is the gravi- 
tational collapse of the core of a massive progenitor star 
exceeding its Chandrasekhar mass. The induced electron 
captures, reducing the degeneracy pressure of the core, 
enable the collapse to proceed until the matter density 
is high enough for nuclear forces to become repulsive. 
This is the case when, roughly, nuclear matter satura- 
tion density has been reached. At this point, a bounce 
occurs, leaving a compact remnant : a proto-neutron star 
(PNS), or possibly a black hole (BH). 

The bounce creates a shock that propagates outwards 
and soon stalls, having lost a lot of energy in photodis- 
sociation of iron group nuclei within the infalling mate- 
rial. The shock is known to stall in the semi-transparent 
regime, where the neutrinos are decoupled from the fluid 
but can still interact. In particular, they deposit en- 
ergy behind the shock via charged-current interactions. 
While the collapse, bounce and prompt shock propaga- 
tion phases seem to be well reproduced by spherically 
symmetric simulations, it is now widely admitted that 
the late phases are deeply multidimensional. 

Over decades much effort has been devoted to numer- 
ical simulations in order to gain physical insight. In- 
creasing computer power and refined models have lead, in 
recent years, to considerable progress in understanding, 
in particular, the supernova explosion mechanism using 
multidimensional models. Indeed, aided by convection 
(e.g. m) and the standing accretion shock instability 
(SASI) [a, Q, several authors reported explosions by the 
neutrino heating mechanism in 2D or 3D (e.g. |7H10l|). 
Although some drawbacks still exist, for instance the re- 
ported explosion energies are low (~ 10 50 erg compared 



with the canonical observed value of ~ 10 51 erg), reviving 
the shock by depositing neutrino energy in the gain layer 
in a multidimensional simulation seems to be a promising 
way to finally make the supernovae explode (e.g. (3,[Tl|). 
However, other mechanisms like the acoustic mechanism 
(e.g. [!,[lH) or the QCD phase transition mechanism [Tp 
cannot be excluded. The MHD mechanism (e.g. |14l - tl7| ) 
could also account for the explosion of rapidly rotating, 
highly magnetized cores and lead to luminous supernovae 
or hypernovae. 

Black hole formation has been intensely investigated 
(see e.g. (l8l - [2l| ). A simulation in which the shock is not 
able to recover positive velocities and to break through 
the infalling material from the progenitor leads to the for- 
mation of a BH, which often swallows the entire progen- 
itor [22j]. This is often referred to as a failed supernova. 
Neutrino emission in failed supernovae simulations stops 
abruptly when the emission region enters the apparent 
horizon, as one could expect. This gives a good crite- 
rion to discriminate between BH formation and neutron 
star formation in an upcoming neutrino signal [23]. Due 
to larger accretion rates, failed supernovae simulations 
are also known to explore higher densities and temper- 
atures than their exploding counterparts, which makes 
them a tool to explore finite temperature equations of 
state (EoSs) at supranuclear density. 

The EoS remains one of the uncertainties in a stellar 
core-collapse numerical simulation. The vast majority of 
simulations employs either the H.Shen et al. EoS (HShcn, 
H) or the Lattimer and Swesty EoS (LS, [25]). Both of 
them assume the same particle content: free nucleons, a- 
particles, one (representative) heavy nucleus, electrons, 
positrons and photons. Concerning the nuclear interac- 
tion, the HShen EoS is based on a relativistic mean field 
model with a Thomas-Fermi approximation for the de- 
scription of the inhomogeneous part, while the LS EoS 
uses a non-relativistic Skyrme-type model with a sim- 
plified momentum-independent nucleon-nucleon interac- 
tion. The latter EoS is further described in Sec. Ill El EOS 
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dependence has been studied since many years with, in 
particular, studies of the effects of incompressibility in LS 
EoS [26j |. or differences between HShen and LS for BH 
formation [l8l - l20j . More recently, it has been shown [27j 
that differences in the nuclear properties can account for 
differences in the resulting hydrodynamics, specifically 
the predicted central density, pressure or compactness of 
the proto-neutron star. Suwa et al [27jj have also shown 
that there is a difference in the gain layer where convec- 
tion is triggered and the simulations behave differently, 
qualitatively and quantitatively. These differences can 
be large enough to govern the presence of an explosion. 

Some other finite temperature EoSs are available, es- 
pecially recent ones (e.g. [2l|, [28j]), focusing on improve- 
ments in the subsaturation regime. In this paper, we 
focus on a different aspect, not included in the standard 
EoSs, namely the presence of additional particles at high 
densities and temperatures. It is natural and widely ac- 
cepted that additional particles should appear in mat- 
ter at densities above roughly nuclear matter saturation 
density [29| . The high temperatures during core-collapse 
are in favor of the population of these additional states. 
Supernovae simulations with pions and hyperons added 
to the HShen EoS were reported in Ishizuka et al. [30| . 
Sumiyoshi et al (3l| and Nakazato et al. [32| . The re- 
sults show in particular that the time from bounce to 
BH collapse is shortened by the presence of pions and/or 
hyperons. The neutrino signal is hardly, or not at all 
changed, except for the fact that the duration of the sig- 
nal is shorter due to earlier collapse to a BH. The reason 
is that pions and hyperons appear only deep inside the 
PNS. Nakazato et al. [32| show in addition that a vari- 
ation of the coupling parameters of hyperons have only 
very little influence on these qualitative results. We ex- 
tend this work to the other commonly used EoS, the LS 
EoS. 

We discuss two extensions of the LS EoS: one including 
a free pion gas and another one including the A-hyperon. 
Unfortunately, very few constraints on the interactions 
of these particles exist, be it from nuclear physics exper- 
iments or astronomical observations. The most stringent 
one at the moment is probably the observation of an al- 
most 2 Mq neutron star [33| , a mass measured with high 
precision by Shapiro delay. Many EoSs including addi- 
tional particles, in particular hyperons, are in contradic- 
tion with this observation due to a too strong softe ning 
of the EoS. For instance, with the HShen EoS + A [34] 
only a maximum neutron star mass of ~ 1.75 Mq is ob- 
tained. With the Ishizuka et al. EoSs [3fl | , the maximum 
masses for their EoSs with different parameterizations of 
hyperons, muons and pions are 1.55,1.63 and 1.65 Mq. 
Our extensions of the LS EoS have been described in 
detail in [35j, where it was shown in particular that it 
is possible, including hyperons in the EoS, to fulfill the 
constraint from the 2 Mq neutron star measurement to- 
gether with the available experimental constraints on hy- 
peron couplings. In Sec. Ill El we shall introduce the EoSs 
models we are employing. 



Our simple neutrino treatment does not allow to ex- 
tract an energy-resolved neutrino signal, we thus only dis- 
cuss some results for integrated luminosities and we will 
concentrate the discussion on the time between bounce 
and black hole collapse as well as on the properties of 
the different PNSs. The EoS with the A-hyperon shows 
a phase transition [36| to hyperonic matter, whose im- 
prints on the collapse shall be discussed. 

This paper is organized as follows. In Sec. |Hl we re- 
view the equations solved by our code and describe the 
employed finite temperature EoSs. In Sec. IIIII we de- 
scribe the newly implemented leakage scheme, and then 
our results are discussed in Sec. IIVI 

In this paper we use a metric signature (—,+,+,+) 
and geometrical units in which c = G = ks = 1- Greek 
indices run from to 3, while Latin indices run from 1 
to 3. We adopt the Einstein summation convention, too. 



II. MODEL AND EQUATIONS 

We perform a series of numerical simulations using the 
CoCoNuT code in spherical symmetry. Note that the 
code can be run in 2D/3D too, but due to the very long 
simulation times, we restrict ourselves to ID simulations. 
Nevertheless, the model and equations presented here 
could be in principle applied to 2D or 3D cases, without 
any loss of gener ality . Since details are already available 
in the literature (37H39| . in this section we only briefly 
describe the equations solved by the code. The newly 
implemented neutrino treatment is described in detail in 
the next section. 



A. General relativistic hydrodynamics 

General relativity (GR) is used through the 3+1 ap- 
proach (see e.g. [40( and references therein), where the 
4-metric g^ v is described by the lapse function a, the 
shift 3- vector (3 l and the spatial 3-metric 7y. The line 
element is then written as 

ds 2 = -a 2 dt 2 + -fij{dx l + P i dt){dx j + /3 j dt) . (1) 

Here we use the isotropic gauge, in which the 3-metric 
is written as 

Hi = - (2) 

with fij denoting the flat-space 3-metric and $ the con- 
formal factor. Although this gauge choice is fully valid in 
spherical symmetry, it is denoted as conformally flat con- 
dition (CFC, see [4l|) in the general case (2D or 3D). This 
approximation misses some aspects of GR (it does not 
contain gravitational waves, it cannot exactly describe a 
Kerr black hole or a rotating fluid configuration), but it 
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has been proven to be a very good a ppr oximation in the 
case of core-collapse simulations [42l |43|. 

Assuming here a perfect fluid with 4-velocity u M and 
baryon number density n.&, we can define a 4-current as 
ribU^. For convenience, CoCoNuT works with the fluid 
density p — m n nb, where m n is the neutron mass, so that 
the 4-current reads 



therefore described as 3-tensors (i.e. using Latin indices) 
tangent to the spacelike hypersurfaces. 

The Eulerian observer sees the fluid with a 3- velocity 
v % = u l I '{ecu ) + f3 l I a. In order to simplify notations, we 
also introduce the velocity v l = v l — ft 1 /a. To write ex- 
plicitely the hydrodynamic equations in flux-conservative 
form [45j . we define a set of "conserved" variables 



J" = pu» . 
The stress-energy tensor is given by 
T^ v = phu"u v + Pg^ 



(3) D = pW, S l = phW 2 v\ £ = E-D = phW 2 



(4) 



with P being the fluid pressure and h — (e + P)/p the 
specific enthalpy; e is here the energy density of the fluid. 
These two quantities obey local conservation equations 



and 



, 



(5) 



where V M denotes the covariant derivative associated to 
the 4-metric g^ v . The electron fraction Y e = n e /rib, 
electron neutrino fraction Y Ue = n l , i! /nb and electron an- 
tineutrino fraction V Ce = n Pii /nb enter separately three 
advection equations, which are included hereafter in a 
conservative-like form. 

No advection equation has been introduced for Y"„ x , 
with v x — v T , Vfj,, u T } standing for muon and tau neu- 
trinos. Before bounce, Y Ux vanishes to a very good ap- 
proximation, because of the absence of charged-current 
production processes due to the lack of heavy charged 
leptons. The only production channels are pair produc- 
tions (see Sec. IIII B 2[) which start to play a role in the 
hot environment close to and after bounce. After bounce, 
we determine Y Um from the integration of the Fermi-Dirac 
distribution with zero chemical potential for the v x , thus 
assuming thermally equilibrated v x . This simplification 
saves some computational time and should not have any 
major influence on our results, since non-electronic neu- 
trinos are much less abundant than electronic ones. Non- 
electronic neutrinos have thus only little effect on the 
overall dynamics compared with the electronic ones (see 
e.g. where the dynamics of a full Boltzmann simu- 
lation is well reproduced by a simplified simulation with 
no v x ). 

The 3+1 formalism slices the four-dimensional space- 
time by three-dimensional spacelike hypersurfaces, and 
it is therefore particularly well adapted to describe hy- 
drodynamics with respect to a Eulerian observer, which 
moves orthogonally to these spacelike hypersurfaces with 
a 4-velocity n M . The 3 + 1 decomposition of the stress- 
energy tensor reads 



T P» _ En^n v + n^S v + S»n u + , 



(6) 



E being the matter energy density, S l the matter mo- 
mentum density, and S lJ the matter stress tensor, all 
measured by the Eulerian observer. Note that the mo- 
mentum density and the matter stress tensors are orthog- 
onal to the Eulerian observer's 4-velocity n M , they can be 



P-D, 
(7) 

where W = au° = l/yjl — ViV 1 is the Lorentz factor. 

From Eqs. ([5]) and (O, adding the three advec- 
tion equations for Y e , Y Ue , Yn. we then derive flux- 
conservative hyperbolic Eqs. [451 ]: 



dt 



dx l 



(8) 



where 



U = [D,Sj,e,DY e ,DY Ue ,DY Ve ] , 

F l = [Dv^SjV* + &)P 1 £v l + Pv\DY e v\DY Ve v\DY^v l ] 

dx l 

da 
dx l 



g and 7 are the determinant of g^ u and 7^,,, respec- 
tively, and Kij is the extrinsic curvature, defined by 
Eq. (TT2"1) . E e , T, Ve , S Pc are the creation-annihilation terms 
that shall be detailed hereafter. These general-relativistic 
hydrodynamic equations are solved using high-resolution 
shock-capturing schemes (see e.g. (4(| and references 
therein). 

The creation-annihilation terms are computed as fol- 
low: 



(9) 



E„. = S 



ec \ ^ee 



Zpl ~ 7^ , (10) 

Y 

Ep C + E ee + Yipi — - — , (11) 



where E ec denotes the electron capture rate, S pc the 
positron capture rate, S ee and E p ; the pair annihila- 
tion and plasmon decay rates, respectively (defined in 
Appendix [Bj . We have introduced here the neutrino leak 
term Y v jt esc that represents the number of neutrinos able 
to reach the neutrinosphere, i.e. able to escape. t esc is 
the time needed for the neutrinos present at a given ra- 
dius r to reach the neutrinosphere and is discussed in 
Sec. IIII B 21 too. 

Note the following modifications compared with previ- 
ous work [381 ]: 
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1 . The neutrino pressure no longer appears in the hy- 
drodynamic source terms, Eq. (|S]). It is instead 
taken into account together with the other neutrino 
source terms, see Eq. (I53"|) . 

2. Two equations for Y Ve and Yp e are added to describe 
neutrino advection. 

3. Three sources, £ e ,Ei, e and S Pc , are added. T, Ue 
and Sp e track creation and annihilation of neutri- 
nos and, via lepton number conservation, S e ac- 
counts for the changes in the electron fraction. 

All these are further detailed in the Sec. Mil and explicit 
expressions are given in App. [3] and App. [B] Simulations 
are performed on a Eulerian grid with spherical geometry. 
Our spherically symmetric runs have 805 (r) points, and 
a logarithmic grid spacing. 

B. Gravitational potential equations in isotropic 
gauge 

Our treatment of the gravitational potential is similar 
to [13, ■ We briefly review here the set of equations 
solved by lorene [47] . a C++ library designed to solve 
3+1 GR equations by means of spectral methods [isjj . 
The CoCoNuT code uses this library for solving gravi- 
tational field equations through the "marriage des mail- 
lages" approach (37} • 

We start by defining the extrinsic curvature 

Kij = + vipi ~ \ fl3Vkpk ) ' ( 12 ) 

where T>j is the covariant derivative associated with 
the flat metric fy. We can then introduce the quantity 

A ij = $ 10 K ij , (13) 
and the momentum constraint can be written as 

VjA ij = 8vr$ 10 S* 1 = 8vr$ 6 .f j Sj = 87r/' iJ S**, (14) 

with S* = $ 6 5j. 
Now, we write as 

A ij = V l X ] + V ] X l - -V k X k f ij , (15) 
3 

where X 1 is a 3-vector (see (39[, too). Combining 
Eqs. (HH) and $T5\) we get 

AX 1 + \v l V. } Xi = torf'Sj, (16) 

o 

which is the first equation solved by the code to obtain 
X % \ we can then deduce A 1 ^ thanks to Eq. (Ti"5"|) and solve 
the Hamiltonian constraint equation 



f.,f. Ahn Aij 

A<t> = -2Trt>- 1 E* -<5>- j IllIjm ^ , (17) 

8 

where E* — Q e E. Knowing the conformal factor, we 
can solve 

A($a) = 2wa$- 1 (E* + 25*) + Qr 7Wi!^! ] 

(18) 

with S* = <f> 6 -yijS ij = $ 6 (ph(W 2 - 1) + 3P) the con- 
formally rescaled trace of the matter stress tensor. It is 
thus possible to get the lapse. Then, by solving 

A/3 4 + -V\V 3 ^) = Vj(2a<S>- e A lJ ), (19) 
3 

we obtain the shift. Note that this approach is not 
bound to the spherically-symmetric case (see [39l | in the 
3D and CFC case). 



C. Apparent horizon 

We use the apparent horizon (AH) finder described in 
details in [49| . In a simulation, a PNS accretes mass and 
looses energy due to neutrinos at the same time. This en- 
ables the PNS to further contract, and consequently the 
density increases. At some point, depending on the EoS, 
the PNS becomes gravitationally unstable and collapses 
to a BH. The BH collapse is unambiguous, as density, 
temperature, pressure, internal energy all increase very 
rapidly. Since the lapse is decreasing from a ~ 1 in 
the Newtonian limit to zero at the center, some authors 
just take a given value a <C 1 of the lapse as a criterion 
for BH formation (see e.g. [H|]). However, using an AH 
finder that detects a marginally trapped surface is a more 
rigorous way to treat BH formation. In addition, the AH 
finder enters a simulation only at the very end, when BH 
collapse begun, triggered by a chosen value for the con- 
formal factor at the center (typically, $(r = 0) > 2), it 
is thus computationally cheap. 



D. Initial models 

We employ two different 40 M Q (zero- age main se- 
quence (ZAMS) mass) progenitors. As a first progenitor 
model, we have chosen a star with low metalicity (10 
solar) from Woosley, Heger, and Weaver [50( (hereafter 
u40). The low metalicity leads to a much higher accre- 
tion rate after bounce compared with a solar metalic- 
ity star. Therefore, even though there still are many 
unknowns concerning the fate of a particular progenitor 
(see e.g. [22|). very massive stars with low metalicity are 
widely accepted as good candidates for BH formation. 
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The second one, from Woosley and Weaver [51J (here- 
after WWs40), has solar metalicity. This progenitor was 
widely used in BH collapse studies (see e.g. [l8l - l2ll |) and 
we use it here mainly for comparison purposes. It is be- 
lieved to collapse to a BH, as the accretion rate is too 
high and the iron core is too massive for the star to ex- 
plode. Let us stress, however, that Ugliano et al [22[ 
report that the 40 M© progenitor with solar metalicity 
from Woosley, Heger, and Weaver does explode and 
leaves a NS rather than a BH. The reason is probably 
that it has a smaller accretion rate than its counterpart 
from [5lj . a difference that may come from the different 
treatment of the mass loss. 



E. Equations of state 

In the central region, i.e. mainly within the hot PNS, 
very high densities (roughly above nuclear saturation) 
and temperatures (several tens of MeV) are reached. Un- 
der these conditions, the particle content of the standard 
nuclear EoS such as the HShen EoS or the LS EoS is 
possibly not sufficient. As already mentioned, the latter 
models the matter as a mixture of one (average) heavy 
nucleus, a-particles, free nucleons, electrons, positrons 
and photons. Additional particles, such as thermal pions, 
hyperons, or nuclear resonances should appear. Even a 
QCD phase transition is conceivable [13]. In the last 
years some work in this direction has begun, extend- 
ing the HShen EoS to include pions or hyperons, see 
e.g. [Mill HI, or quarks, see [ij]. 

Here we employ EoSs which are extensions of the other 
widely used nuclear EoS by Lattimer and Swesty [25| . 
Electrons and positrons are treated as non-interacting 
relativistic gas in pair equilibrium, neglecting electron- 
screening effects; photons are treated as an ideal ultra- 
rclativistic gas. Equilibrium with respect to strong and 
electromagnetic interactions is assumed, while j3 equilib- 
rium is not requested, which is consistent with expecta- 
tions during core-collapse events. Nuclear interaction is 
treated using a liquid drop model and the transition be- 
tween inhomogeneous and homogeneous nuclear matter 
is described via a Maxwell construction. Details can be 
found in the original paper 25] . 

There are three different sets of parameter values 
available, resulting in three different values of the nu- 
clear incompressibility, K = 180, 220, and 375 MeV. In 
this work, we restrict ourselves to the value K — 220 
MeV, which gives an EoS compatible with both nu- 
clear data, which suggest a value around K — 240MeV 
(see e.g. [52]). and the recent high mass neutron star of 
M = 1.97 ± 0.04 M Q [H]. The maximum mass of a cold 
neutron star is 2.06 M© (2.40 M© baryonic mass) with 
K = 220 MeV. The EoS with K = 180 MeV, on the 
lower end of possible X-values from nuclear data, fails 
to reproduce the neutron star mass constraint giving a 
maximum cold neutron star mass of only 1.84 M© and 
K = 375 MeV lies well above the allowed region from 



nuclear data. For a thorough comparison of the three 
different LS EoSs in the context of core collapse simula- 
tions, see [27l ] . 

In [35{ different possible extensions of the LS EoS are 
discussed, considering hyperonic degrees of freedom as 
well as pions and muons. Within this work we limit the 
discussion to two rather simple cases: one including pi- 
ons and the other one including the A-hyperon. For the 
former one, pions have been added upon the LS EoS as 
a fr ee g as. For a critical discussion of this approximation 
see [35| . 

The second one, including the A-hyperon can be seen 
as analogue to the HShen + A EoS from [34]. The mo- 
tivation is that the A represents, together with the E~- 
hyperon, probably the most important hyperonic degree 
of freedom in hot dense supernova matter. Thus, includ- 
ing the A allows for discussing general features of the 
effects coming from the hyperonic degrees of freedom, 
without the necessity of resolving the complicated parti- 
cle composition in the presence of many different hyper- 
ons. Indeed, most of them have very low abundances, see 
e.g. (35[. The AA and the AA^ interactions are taken from 
the model by Balberg and Gal [53[ with the parameter- 
ization 220BG from 35J. This model has the advantage 
of matching well with the LS EoS: in the region where 
no A- hyperons are present, the interaction is exactly the 
same as in the LS model. 

Within the Balberg and Gal model, the A-hyperons ap- 
pear through a first order phase transition, see |36l. l54l] . A 
Gibbs construction is employed to describe matter in the 
phase coexistence region, see [36| for details. The effect 
of the criticality of the phase transition on the simulation 
shall be the subject of forthcoming work. 

The maximum gravitational masses of cold neutron 
stars are M = 1.95 M Q for the LS220+7T EoS (2.27 M Q 
baryonic mass) and M = 1.91 M for the LS220+A EoS 
(2.22 Mq baryonic mass) with parameterization 220BG. 
The former value lies in the ler-range for the mass con- 
straint from PSR J 1614-2230, the latter is only very 
slightly below. Properties for cold neutron stars for the 
different EoSs used within this work are summarized in 
Table. D The LS model with K = 180 MeV is added for 
comparison. The circumferential radii at the canonical 
value of M — 1.4 M Q vary only very little for the four 
EoSs, whereas the radii at maximum mass are lower for 
the EoSs with additional particles, as expected from the 
softer character of these EoSs. 



III. LEAKAGE SCHEME 

A. Introduction 

Neutrino transport is one of the most challenging as- 
pects of modern supernovae simulations. Because the 
shock stalls in the semi-transparent regime, an accu- 
rate neutrino treatment should rely, in principle, on solv- 
ing the Boltzmann equation. Several authors report on 
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TABLE I: Properties of cold neutron stars for the different employed EoSs. The LS180 EoS is added for comparison. The given 
radii are circumferential ones. 



Name 


LS180 LS220 LS2204 


A LS220 + ix 


Maximum gravitational mass [Mq] 


1.84 


2.06 


1.91 


1.95 


Maximum baryonic mass [Mq] 


2.12 


2.40 


2.22 


2.27 


Radius at maximum mass [km] 


10.13 


10.67 


9.28 


9.94 


Radius at M a = 1.4M Q [km] 


12.19 


12.71 


12.41 


12.03 


Central density at maximum mass [fm -3 ] 


1.26 


1.09 


1.47 


1.26 



supernovae simulations with neutrino transport, cither 
Newtonian (e.g. [H, [5(|), or general relativistic (e.g. 
[57l |58| ) , but this remains computationally challenging : 
a simulation with 3D GR hydrodynamics and a (6D) 
Boltzmann solver is not done yet (see also 23] ) . The use 
of simplified models is therefore fairly common. They 
are used to tackle 3D hydrodynamics simulations [59j . 
or for parametric studies of black hole collapse, which is 
known not to require a fully detailed neutrino transport 
[2fl [22| . In this section, we provide details concerning the 
implementation of our leakage scheme, mostly inspired by 

The main idea of the leakage scheme is to treat the 
neutrinos as a fluid component inside the neutrinosphere. 
Above the latter, neutrinos are considered as free stream- 
ing, and only the energy taken away from the fluid needs 
to be considered. This is achieved by adding terms to the 
sources of hydrodynamic equations, Eqs. ([8]). These addi- 
tional source terms are explicitly given in Appendix [XJ It 
is a grey scheme, that defines and works with a mean en- 
ergy for every neutrino species. The leakage scheme is, by 
construction, not sufficient to treat the semi-transparent 
regime nor to revive the shock. On the other hand, it is 
a simplified model that comes from physical arguments 
and enables us to run core-collapse simulations within a 
reasonable amount of time 1 . 



B. Implementation 

1 . Neutrinosphere 

The neutrinosphere is the limit where neutrinos de- 
couple from the fluid. This is characterized by an optical 
depth t of the order of unity. We choose to define it 
as the region where r = 2/3, consistently with the litera- 
ture (see e.g. [IHIgij]). In order to find the exact position 
of the neutrinosphere during a simulation, at each time 
step the opacity is evaluated, using the mean energy of 
the neutrinos, Eq. (l37l) . 



Following [TH 1601 - 1621 ] , the neutrino reactions we con- 
sider for calculating the opacity are the following (see 
also App. |B]for the details of the implementation) : 

1. Elastic scattering off a nucleon, 

Vi + N -^Ui + N , (20) 

where N — n,p, and Vi represents one of the three 
neutrino species implemented : electron neutrino 
v e , electron antineutrino v e or other species v x . 

2. Elastic scattering off a nucleus, 

Vi + (A,Z)->Vi + (A,Z), (21) 
where (A, Z) is the mean nucleus. 

3. Absorption of a v e by a neutron 

v e +n-tp + e~~ (22) 
and absorption of a v e by a proton 

v e +p^n + e+, (23) 

are taken into account, when no /^-equilibrium is 
assumed. 

4. Elastic scattering off a A 

Vi + A Vi + A . (24) 

We implement the coherent scattering off A- 
hyperons and take it into account in the opacity 
calculation, because the hyperon fraction can be- 
come quite large at the end of a simulation (see 
Sec. |TV] for the hyperon fraction and App. [B] for 
the implementation). We study the impact of this 
newly implemented reaction on our simulations in 
Sec. HVC3l 

The optical depth is obtained from the total opacity, k, 
by integrating over r, 



1 e.g. fiducial model lsu40 takes 152 minutes for a single-processor 
run. The total physical time simulated in the code is 871ms. 



k dr 



(25) 
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Note that this expression implicitly assumes that neutri- 
nos only move along radial rays. In our code, the upper 
bound of the integral is taken to be the last cell, where 
we start the integration, then we move inwards. The first 
cell reaching r > 2/3 is taken to be the neutrinosphere. 

As a refinement of this basic leakage scheme, an ef- 
fective neutrino chemical potential Hv,eff is introduced 
in the calculation of the different cross sections follow- 
ing [13, 

fiv t eff = A^,eg(l - exp(-r)) , (26) 

with fj, u ,eg being the neutrino chemical potential in /3- 
equilibrium given by the EoS as Hv,eq — Me + Mp — Mn- 
For matter completely transparent to neutrinos with a 
vanishing optical depth, the effective chemical potential 
goes to zero. In the other extreme case of a very large 
optical depth, corresponding to matter opaque to neu- 
trinos, the equilibrium value is reached. Thus, it allows 
partly to correct for deviations from equilibrium in the 
semi-transparent regime, although the neutrino distribu- 
tion function is still represented by a Fermi-Dirac distri- 
bution. 

From a practical point of view, Eq. (|26l) usually needs 
very few iterations to converge (no more than 3), given 
the neutrino effective chemical potential from the pre- 
vious time step as a first guess. At the beginning of a 
simulation, (J, v ,eff is taken to vanish everywhere. 

2. Neutrino creation and advection 

Neutrinos are created by electron capture on free pro- 
tons (see App. [B]for the details of the implementation) 

p + e~^v e +n, (27) 

electron capture on nuclei 

(A,Z)+e- -> (A,Z-l) + u e , (28) 

positron capture on free neutrons 

n + e + ^P e +p, (29) 

electron-positron pair annihilation 

eT + e+ ->• Vi + P i: (30) 

and plasmon decay 

7 ->• i>i + P^ (31) 

Since they appear only in the very dense regions close to 
the center, where /3-equilibrium is assumed, no charged- 
current reactions on A-hyperons are considered. 

Let us stress several approximations applied in the 
treatment of electron and positron captures. We can jus- 
tify these approximations by the fact that the leakage 
by itself is a rather crude approximation, and going into 



further details without a better neutrino treatment may 
not be relevant or possible. 

Following the rates given in Bruenn [gij , electron cap- 
ture on nuclei is cut off due to the assumption of shell clo- 
sure for a neutron number N > 40. It has been shown [65| 
that temperature effects smear out the shell structure and 
that the capture rate on nuclei with N > 40 is therefore 
nonzero. 

Electron and positron capture are only relevant out 
of /3-equilibrium, since in equilibrium the creation and 
absorption terms cancel each other. Within the leakage 
scheme, however, the equilibration cannot be described 
as the balance between both would require the knowledge 
of the distribution function. We therefore take only the 
creation part into account below a given critical density 
where /3-equilibrium is assumed. We choose this density 
to be 1.17 x 10 12 g/cm 3 (see Sec- lmD]) . which is m agree- 
ment with the general finding that /3-equilibrium sets in 
at a density of the order p ~ 10 12 g/cm 3 , see Sec. IIII B~3l 
too. 

In the region behind the prompt shock, matter is 
strongly neutronized, giving very low values of Y e . This 
effect is in general overestimated by leakage schemes (see 
Fig. [2 or e.g. (62[), among others due to the fact that the 
absorption term is absent, which is slowing down electron 
capture in this region close to /3-equilibrium. In order to 
avoid the leakage scheme giving (unrealistic) values of Y e 
below the lowest one available in the EoS, electron cap- 
ture is blocked below Y e = 0.045. This allows to follow a 
simulation during and after the prompt shock propaga- 
tion 

The reaction rates for neutrino creation and absorption 
enter the sources of the new advection equations related 
to electron and neutrino number in Eqs. (f^- fTTj) . They 
are detailed in App. [Bj We give here the expression for 
the neutrino escape time t esc , defined as in [60(: 

t esc = ai(R v - r)r . (32) 

R v is here the radius of the neutrinosphere. Ruffert et 
al. [60| have set the free parameter a\ to a\ = 3, in order 
to reproduce well data from transport calculations. We 
keep this value in our implementation (see Sec. IIII D[) . 

Since we determine the position of the neutrinosphere 
at each time step, some neutrinos may be no longer 
trapped at a given time step. We separately account 
for these freed neutrinos due to the inward movement of 
the neutrinosphere. The treatment of energy losses due 
to escaping neutrinos will be presented in Sec. IIII CI 

Finally, note that Eq. (|3"2")l as well as Eq. (|2"5]l do not 
have GR corrections. Neutrinos follow a straight line 
(on the grid) while in principle they should follow a null 
geodesic. This introduces a small numerical error, com- 
pared with the overall approximations induced by the 
usage of a leakage scheme. In particular, the ad hoc fac- 
tor of a\ = 3 in Eq. (j32|) is already larger than the GR 
corrections in Eqs. (|32l25p during our simulations. The 
biggest underestimation of opacities occurs during BH 
formation, because the neutrinospheres move closer to 
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the center (or to the newly formed apparent horizon). 
But we note that we are able to reproduce the collapse 
of the PNS, which is essentially a free fall. 



3. Deviation from thermal equilibrium 

During a core collapse event, neutrinos are in general 
not in thermal equilibrium, which results in a distribu- 
tion function different from the Fermi-Dirac one, when 
using a realistic neutrino transport scheme. Only in some 
very small regions it can be approximated by the latter. 
The simplicity of the leakage does not allow for determin- 
ing the correct distribution function, implying that the 
neutrino number given by integrating the Fermi-Dirac 
distribution function is not the same as that given by 
advection and creation- absorption. Within some leakage 
schemes the neutrino number is chosen to be given by 
integration of the Fermi-Dirac distribution function, as- 
suming thus thermal equilibrium at all times. We have 
chosen to compute instead the neutrino fractions from 
source terms and advection equations. Therefore, the 
current neutrino number is not necessarily the same as 
the value at thermal equilibrium, which is taken only as 
a maximal value, similarly to [63[. 

This choice allows to perform simulations through the 
collapse and bounce phases, where otherwise the neu- 
trino number would be strongly overestimated. To illus- 
trate this, in Fig. Q] we display the F e -pronle at bounce, a 
few milliseconds before the launch of the shock. The re- 
sult obtained with our leakage scheme is compared with 
that employing the so-called Liebendorfer deleptoniza- 
tion (see [66|). The latter is a parameterization of full 
Boltzmann simulations and reproduces thus well a real- 
istic profile. We can see that the result of the leakage 
is in good agreement with the predicted deleptonization. 
The deviation is at most around 5%, excepted at the 
very center where electron captures in the Liebendorfer 
scheme are slowed down and finally stopped due to the 
approaching and onset of /3-equilibrium. Y e becomes con- 
stant at a value around 0.28. The leakage scheme, which 
neglects the absorption term for electron capture can- 
not reproduce this trend. Y e continues to decrease un- 
til the critical density for the onset of /3-equilibrium is 
reached. Since this density value is adjusted to repro- 
duce the Ye-profile as well as possible during neutroniza- 
tion after bounce (see Fig. [5] and discussion below), Y e 
at the center at bounce is a little low (Y e w 0.22), but 
this should not have a significant impact on the overall 
dynamics (see [63], too). 

Fig. [2] shows a comparison of the F e -profile in the leak- 
age scheme with that from the code Agile/Boltztran 
(from Liebendorfer et al. [HI) with a Boltzmann neutrino 
solver at 10 ms after bounce. This corresponds to the pe- 
riod where matter behind the propagating prompt shock 
is strongly neutronized. Note that, in order to have the 
same setting as [68j , we show here a spherically symmet- 
ric simulation with a 15 M Q progenitor from Woosley and 
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FIG. 1: Ye as a function of radius at bounce, using the leak- 
age scheme (plain, red) compared to the Liebendorfer delep- 
tonization prescription (dashed, green) 
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FIG. 2: Yo-profile during neutronization at t — 10 ms post- 
bounce using the leakage scheme (plain, red) compared with 
Agile/Boltztran (dashed, green). 



Weaver [51| with solar metalicity. Moreover, in contrast 
to all other runs we have performed, the EoS used here is 
the LS EoS with an incompressibility of K = 180 MeV, 
to match the setup of the run with Agile/Boltztran 
from j68| . 

We observe that, as mentioned above, one of the draw- 
backs of the leakage scheme is to overestimate the neu- 
tronization. The main reason is that the equilibration of 
electron capture is not correctly described. Otherwise, 
the leakage scheme is in qualitative agreement with Ag- 
ile/Boltztran. Note that the irregularities in Y e are 
due to the different simplifications of the leakage scheme. 
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C. Treatment of energy losses 

In this section, we describe how the energy lost by the 
fluid, because of neutrino emission, is taken into account. 
No heating is implemented, as the leakage is not a trans- 
port scheme, and therefore no self-consistent derivation 
of heating sources exist. Moreover, heating does not play 
a major role in our simulations, given that we expect a 
collapse to a BH (see discussion in Sec. Ill D|) . 

Neutrinos enter the fluid momentum and energy equa- 
tions as a source term in the equation VpT a @ = 0, 

V p T af3 =q a , (33) 
where we define the source 4- vector q a as 



Qlf = 



( Qe \ 

Q T (M) 

Q 



\ 



(M) 
(M) ) 



(34) 



If the full neutrino distribution function were known as 
the solution of the Boltzmann equation, the latter could 
be employed to build a consistent neutrino stress-energy 
tensor from which the hydrodynamic sources q a could 
be derived relying on the condition that the total stress- 
energy tensor is divergence-free. This has been done in 
GR by several authors (see e.g. [57], [H, [6(|), incorporat- 
ing the energy losses in a coherent way. This is, however, 
not possible in a leakage scheme, where the energy losses 
can only be implemented in an approximate way. We fol- 
low the choice of [62[ for the vector q a , which is detailed 
below. The quantities in Eq. (f3"4"|) are defined in the fluid 
rest frame (or Lagrangian frame, hence the subscript LF), 
we further need to transform them to the Eulcrian frame 
where the hydrodynamic equations, Eqs. ([8]), are solved. 
This transformations are detailed in Appendix [Al where 
the explicit expressions for the source terms used in the 
code are given as well. 



1. Free streaming regime 

In the free streaming regime all neutrinos that are pro- 
duced leave the simulation, and take away their energy 
from the fluid. Knowing the neutrino creation rates, Ri, 
from each considered reaction, and the mean neutrino 
energy in each reaction, (e^), (following again [I(|), we 
can write the total rate of energy loss by the fluid in the 
free streaming regime as follows: 



Qe — S ec (e ec ) + Yi pc (e pc } + S ee (e ee ) + Y]p[(€pi} + Q freed ■ 

We have added here to the energy loss by neutrino cre- 
ation in electron/positron capture, pair annihilation and 



plasmon decay, the energy loss rate, Q freed, due to the 
freed neutrinos (if any) when the neutrinosphere moves 
inwards. This term is new and comes from our imple- 
mentation of advection equations for neutrino fractions. 
See App. [B] for further details on the implementation. 
Note that no momentum is removed from the fluid in the 
free streaming regime. 



2. Trapped regime 

In the trapped regime only those neutrinos reaching 
the neutrinosphere take away energy. The total energy 
per unit time lost by the fluid in the trapped regime thus 
reads 



!E = -a 2 {e v ) 



(36) 



with the mean energy per particle of escaping neutrinos, 
(e„). In our grey scheme, this mean energy per particle 
can be evaluated as follows, 



Ef PD (v)d 3 p _ t f 3 ( v ) 



J °° f FD (v)d 3 P F 2 {r,y 



(37) 



where fp d is the Fermi-Dirac distribution function and 
F]~ is the Fermi integral of order k, 



F k {v) 



Ax 



1 + cxp(a; — rj) 



(38) 



T is the temperature and rj is the degeneracy param- 
eter (r\ v — fi u /T). Note that since (e v ) depends on 7y„, 
its value is sensitive to the use of the effective chemical 
potential introduced in Eq. (|26l) . 

The energy loss rate, Eq. (|3"6"]l . is calculated separately 
for each neutrino species considered (v e , D e , v x ), each one 
having its own mean energy and its own fraction. 

In order to correct partly for the fact that the leakage 
scheme describes the energy loss only in an approximate 
way, we choose to vary a 2 as a free parameter of order 
unity. It is adjusted to have a mass of the PNS at BH 
formation comparable with the literature [ljl [2(| • In this 
paper, we fix a 2 to the value a 2 — 1.1, see Sec. IIII Dl 

In the trapped part, the momentum removed from the 
fluid is approximated as the gradient of the neutrino pres- 
sure, P v [58l l62l|. To compute it, neutrinos are assumed 
to be an ideal Fermi gas of massless particles up to the 
neutrinosphere. The neutrino pressure then reads 



AttT 4 
3(hc) : 



(39) 



The term proportional to F 3 (r) u ) in Eq. (131)|) takes into 
account the v e pressure, while the one proportional to 
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Fs(— r\v) takes into account the v e pressure. The pressure 
of v x is negligible and therefore not considered. 

The momentum taken away from the fluid is then ob- 
tained by computing the gradient of P v via [38l . l62j | : 

BP 

Q(M ) i = -^l (40) 

This expression shall be used to derive the neutrino 
energy-momentum source terms for the hydrodynamic 
equations in Appendix [A"l 

Note that the components QiM)i are exactly zero in 
the free streaming regime. 

D. Summary of parameter values 

Our leakage scheme has three different parameters. 

• The /3— equilibrium density is taken to be 1.17 x 
10 12 g.cm -3 throughout this paper. Different au- 
thors report j3— equilibrium to set in between 6 x 
10 11 and 2 x 10 12 g.cm -3 (Hill. Our value is in 
agreement with this range. 

• oi, defined by Eq. (j3"2")) . has the numerical value 
ai = 3. We keep the value reported in [60j. 

• d2, defined by Eq. (|36[) . has the numerical value 
<i2 = 1.1. Since this parameter corrects for the in- 
completeness of the leakage scheme, it has to be 
close to unity. It is adjusted to have a mass of the 
PNS at BH formation comparable with the litera- 
ture [H,[20|. 

The quantitative values found in our simulations are 
quite sensitive to these parameters. Namely, a difference 
of ±10% in the parameter a 2 leads to a difference in 
the PNS maximum mass of ~ 2%, which can lead to 
significantly more time spent in the accreting phase. A 
difference of ±10% in the j3— equilibrium density does 
not change the AH detection time, but the maximum 
density at bounce changes by about ±10%. These tests 
were done with the lsu40 model (see Sec. II V A|) . 

We note that the overall behavior described in the re- 
sults (Sec. IIV[) is robust. We keep in mind that our neu- 
trino treatment is simplified and that a complete Boltz- 
mann solver may give different quantitative values. 

IV. RESULTS 

We use either progenitor u40 or WWs40 (see Sec. lIIDD . 
and we employ the LS2 20 Eo S, the LS220±?r EoS or the 
LS220+A EoS (see Sec. Ill Eft . We name our simulations 
with the name of the EoS (Is, tt, A), followed by that of 
the progenitor (u40 or WW40). Names and results are 
summarized in Tab. HI] Simulations are stopped after the 
AH is formed. 



A. Fiducial model 

We start by describing the simulation lsu40, using the 
LS220 EoS with no extra particles and the u40 progenitor 
from [50j. This simulation shall serve as a reference run 
in order to compare with simulations using other EoSs 
and/or the other progenitor. 

Both the collapse and bounce phases are followed us- 
ing the leakage scheme, which is possible because we 
solve advection equations for Y Ve and Yp e as described 
in Sec. MI B 31 The bounce is detected at the moment of 
the shock formation, but, as there is no switch between 
the different schemes before and after bounce, the exact 
moment of bounce detection has very little influence on 
the dynamics (see [7l[ for an example of an implementa- 
tion using a switch at bounce). 

After the bounce, the shock is launched and it reaches 
a distance of about 98 km from the center before it stalls. 
Then, the cooling of the material makes the shock slowly 
recede. At the end of the simulation, before the BH col- 
lapse is triggered, the shock remains about 57 km away 
from the center and its velocity has increased to 0.35c. 
The PNS radius, defined somewhat arbitrarily at the 
density of 10 11 g.cm -3 following [72J, recedes, too: At 
bounce, it is about 75 km, then it shrinks slowly to about 
45 km at the onset of the BH collapse. 

As depicted in Fig. [31 the PNS collapses at 415 ms 
post-bounce, within about 2 ms. The central density 
increases from 4.6 x 10 14 g.cm -3 to 1.5 x 10 15 g.cm -3 and 
the temperature raises from 4.3 x 10 11 K to 1.3 x 10 12 K. 
At 416.8 ms post bounce, the apparent horizon is found 
and the black hole is formed. 

The enclosed mass at the onset of BH collapse, mea- 
sured by evaluating the baryonic mass up to the position 
of the shock, is 2.55 M Q (all masses given in this section 
are baryonic masses). It is a measure (slightly overes- 
timated, because the shock is not exactly at the PNS 
radius) of the PNS maximum mass, above which it col- 
lapses to a BH. This result is consistent with [20j within 
3% (their reported PNS mass is 2.469 M ). 

B. Model with pions and the u40 progenitor 

Adding pions to the EoS slightly softens it, so that we 
expect the PNS maximum mass and the PNS radius to 
be slightly smaller. 

During the entire simulation, 7r - are the most abun- 
dant pions. This is to be expected since negatively 
charged particles are favored within the neutron rich hot 
and dense matter encountered during core-collapse. At 
bounce the maximum of is Y n - = 0.0032. The influ- 
ence of pions at this moment becomes non-negligible. 

In Fig. [31 we can see a consequence of the presence 
of pions on the collapse phase, namely, the central den- 
sity at bounce is slightly lower in model 7ru40 than in 
model lsu40. Indeed, collapse and deleptonization are 
slightly different, and it results in an homologous core 
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FIG. 3: Central density as a function of time post-bounce, for models lsu40 and 7ru40. The time corresponding to the first 
detection of the AH is marked by a vertical line. 



and a nascent PNS larger in model 7ru40. More pre- 
cisely, the central temperature is ~ 2% higher at bounce 
and the difference stays relatively constant during the 
first few tens of milliseconds, and this results in a ~ 3% 
lower central density for model 7ru40 compared to model 
lsu40. 

As can be seen in Fig. [31 the PNS for ir u40 is slightly 
more compressible and the central density increases more 
rapidly. At ~ 200 ms post-bounce, there is a crossing 
point where the central densities of both models are the 
same, then the model with pions predicts a higher density 
until the collapse to the BH. The addition of pions only 
slightly changes the behavior after bounce. 

Fig.|4]shows the pion fraction just before the collapse to 
the BH at t = 356 ms post-bounce. As mentioned above, 
we observe that 7r~ are the most abundant, followed by 
7T° and then 7r + . 

We note that 7r~ is the only species with a fraction 
greater than 2%, and all species are present only in the 
very center of the PNS, up to ~ 30 km. The shape of 
the profile, i.e. the fact that near the center the pion 
fractions are smaller, can be explained by the lower tem- 
perature near the center, too. Pions are mainly thermally 
excited, so that their fractions naturally follow closely the 
temperature profile. The pion fractions are of the same 
order as those reported by [32j using an extended HShen 
EoS. The most abundant species, n~ , has a maximum 
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FIG. 4: Pion fractions as functions of radius at the onset of 
BH collapse (t = 356 ms post-bounce, model 7ru40). The 
shape of the profile, in particular the decrease in the central 
region, is directly correlated with the temperature profile. 

of Y w - ~ 0.13 at the onset of the BH collapse. During 
a simulation, all pion fractions increase, from to their 
maximum value at the BH collapse, following the increase 
of the temperature and the density inside the PNS. 
We find that the apparent horizon is detected at 
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366.3 ms post bounce for a PNS maximum mass of 
2.49 M Q with a PNS radius of 43 km. We note that even 
though pion fractions are not very large (Y^- is the only 
one over 0.10, and only in the last tens of milliseconds 
before BH collapse), this small change in the EoS results 
in non-negligible changes in the dynamics, especially the 
time before BH collapse. 



C. Model with A-hyperons and the u40 progenitor 

1. Contraction of the PNS 

Adding A-hyperons to the EoS results in a more com- 
pact PNS, and the effect is much more pronounced than 
with pions. 

Indeed, in model Au40, the onset of collapse to BH oc- 
curs at 122.8 ms post-bounce and the AH is detected at 
123.9 ms post-bounce, for a maximum mass of the PNS 
of only 2.00 M Q . Note that this value is significantly 
lower than its cold EoS counterpart (LS220+A allows 
for cold neutron stars with a maximum baryonic mass of 
2.22 M Q ). This difference might be attributed to dynam- 
ical effects. In addition, the hot proto-neutron star is not 
entirely in /3-equilibrium, such that the Y e profile is given 
dynamically. Since the EoS depends considerably on Y e , 
this might explain the lower supported mass, too. 

Fig. [5] shows that the central density increases much 
more rapidly than in the fiducial model lsu40. The shock 
propagates only up to 36 km away from the center before 
receding. The newly formed PNS extends to 31 km a 
few milliseconds after bounce and contracts to a radius 
of 17 km very rapidly. Contraction continues slowly all 
along the simulation up to the phase transition (described 
in the next subsection), where the PNS radius suddenly 
decreases from 15 km to 13 km. It further contracts to 
about 10 km at the onset of BH c ollap se. This radius is 
very small for a hot PNS (see e.g. [73| for reported PNS 
radii with another progenitor and the LS180 EoS. Radii 
typically lie between 20 and 70 km). 



2. Phase transition 

The LS220+A EoS contains a first order phase transi- 
tion to hyperonic matter, as discussed in Sec. Ill El which 
we expect to find in a core-collapse simulation. Indeed, 
Fig. [6] is a detailed view of Fig. [5l where a phase transi- 
tion occurs. The density raises from 9.4 x 10 14 g.cm~ 3 
to 11.2 x 10 14 g.cm -3 within less than 2 ms (around 68- 
70 ms post-bounce). 

This phase transition leads to a sudden contraction 
of the PNS on a dynamical time scale, as discussed in 
Sec. IIV C II Further contraction of the PNS is very sim- 
ilar to the contraction of the iron core during the initial 
collapse. 

Fig. [7] shows the A-hyperon fraction YX as a function 
of radius during the phase transition. At finite tempera- 



ture, A-hyperons already appear at bounce, and Y\ keeps 
increasing to reach ~ 0.17 at the center at 68 ms post- 
bounce. Then, the phase transition occurs, and we can 
see in Fig. [7] that Y\ increases from ~ 0.17 to ~ 0.22 at 
the center. This happens within less than 2 ms. After 
the phase transition, Y\ oscillates with the density and 
increases again at BH collapse. 

Fig. [5] shows Y\ at the onset of BH collapse. We see 
that Ya is quite large and can go up to 0.41 at the very 
center of the PNS. Note that the phase transition is a 
consequence of the very large accretion rate of the pro- 
genitor, and other progenitors can give a PNS that col- 
lapses directly to BH before reaching the region of phase 
transition (see Sec. lIVDj) . 

Moreover, Fig. [5] shows oscillations of the PNS excited 
by the collapse at phase transition. These could be fun- 
damental mode radial oscillations of the PNS (see [74{ for 
a study in cold neutron stars, with different EoSs). It is 
not possible to directly compare the oscillation frequency 
of our simulations, which we estimate to be ~ 800 Hz, 



to the data presented in [74|, because i) the authors dis- 
cuss only cold neutron stars, ii) our system accretes mat- 
ter continuously and hi) we use none of the EoSs listed 
in their tables. Nevertheless, the order of magnitude is 
correct and given the uncertainties on our estimations of 
mass and radius of our PNS, we might interpret these os- 
cillations to be the fundamental radial mode of the PNS. 
In addition, we performed several simulations with in- 
creased resolution to rule out a numerical artifact, and 
we found these oscillations to be robust. 



3. Influence of the A-scattering 

To test the influence of our newly implemented isoen- 
ergetic scattering off A (see Sec. IIII BTI and App.[B]), an- 
other simulation is conducted using the same parameters 
and the same EoS, but switching off the A-scattering, 
and we compare it to the simulation with A-scattering. 

Differences are negligible up to the phase transition. 
Then, without scattering off A, the phase transition oc- 
curs 6 ms later, and the PNS maximum mass becomes 
1.96 M Q . This is a reduction of 0.04 M Q compared 
with the case including A-scattering. This reduction of 
the maximum mass suggests that cooling is more impor- 
tant when we reduce the opacity by neglecting scatter- 
ing off A-hyperons, as expected. However, the reduction 
is within the error bars of the code, so this conclusion 
should be taken with some care. We note that although 
A-hyperons are abundant, they appear in a medium that 
is already optically thick. This could explain why the 
A-scattering is not very important. 



4- Neutrino luminosity 

As a complementary tool, we compute the total lumi- 
nosity L v . In spherical symmetry, we approximate it by 
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FIG. 5: Central density as a function of time post-bounce, for model Au40. The time corresponding to the first detection of 
the AH is marked by a vertical line. 



TABLE II: Properties at BH collapse. 



Name 


lsu40 


7TU40 


Au40 


lsWW40 


ttWW40 


AWW40 


Initial model 


u40 


u40 


u40 


WWs40 


WWs40 


WWs40 


EoS 


LS220 


LS220+7T 


LS220+A 


LS220 


LS220+7T 


LS220+A 


AH (ms post-bounce) 


416.8 


366.3 


123.9 


1025.8 


607.2 


274.5 


PNS mass (M ) 


2.55 


2.49 


2.00 


2.73 


2.40 


2.00 


PNS radius (km) 


45 


43 


10 


51 


38 


18 



4-7T 



W 1 



a 



-9 dr. (41) 



Qe being denned in the fluid rest frame, we have to 
transform it to the Eulerian frame with a Lorentz boost 
(with a Eulerian velocity v r /& 2 , see Appendix [A"f . and 
then to the coordinate frame, which is done by multiply- 
ing by the terms in brackets. The generalization to 3D 
would involve a Lorentz boost in an arbitrary direction 
and the complete shift vector. 

The neutrino luminosity can be found by integrating 
the energy rate on an invariant volume element at con- 
stant time, which justifies the y/—g term. This formula 
agrees with [62|,[75j]. Note that the integration at constant 



time is an approximation, since a time delay between the 
different emitting regions should exist. But because the 
overall emitting region is narrow (~ 100 km at most), 
this approximation should be well justified. 

As already reported by different authors (e.g., [3"3 , l75j ). 
the total neutrino luminosity peaks at the neutrino burst, 
a few milliseconds after the bounce, and decreases very 
rapidly afterwards. We observe the same behavior (see 
also Fig. [17)1 and the discussion for the WWs40 progenitor 
in Sec. HVPp , and the luminosity for models lsu40 and 
7ru40 are very similar. Only Au40 shows a different trend. 

The neutrino burst peak luminosity is higher, and be- 
cause the PNS shrinks very rapidly, high neutrino emis- 
sion after bounce is more sustained with the model Au40. 
As a result, the luminosity decrease after the neutrino 
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FIG. 6: Central density as a function of time post-bounce, 
model Au40, zoomed on the phase transition. 
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FIG. 7: Ya evolution during phase transition (model Au40). 
Snapshots are taken at the onset (t = 68.7 ms post-bounce), 
half-way through (t = 69.5 ms post-bounce) and at the end 
of the phase transition (t = 70.2 ms post-bounce). 
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FIG. 8: YA-profile (model Au40) at the onset of BH collapse 
(t — 121 ms post-bounce). 



burst is slower in this model. 

Moreover, after the phase transition the luminosity fol- 
lows the PNS behavior (illustrated by the central density 
plot in Fig. [5j, and oscillates, which makes the PNS loose 
a lot of energy. The differences in luminosity between 
LS220 EoS and AEoS have to be checked with a better 
neutrino treatment to infer the detectability of the LS220 
+ A EoS peculiarities. In addition, the detection of PNS 
modes could be linked to a physical process different from 
the phase transition, and it would be hard to disentangle 
on a neutrino light curve. 

Finally, at BH formation, the emitting region enters 
the horizon and the luminosity drops down abruptly. 



D. Results with a WWs40 progenitor 

Within this section we discuss simulations performed 
using a 40 M Q , solar metalicity progenitor from [5lj ]. 
varying between EoSs LS220, LS220+7T and LS220+A. 
This progen itor has often been used to study BH col- 
lapse [18| - |21| . as its iron core is reported to be very large. 
The main difference with the u40 progenitor from [5(| is 
that the accretion rate is significantly higher in the latter. 

In Fig. [9j we observe the same trend as for progeni- 
tor u40. With the model 7rWW40, the bounce occurs at 
lower density than with the model lsWW40, and after 
275 ms the curves cross, due to the fact that LS220+7T 
EoS allows the PNS to contract faster. For model 
lsWWs40 we find a maximum PNS mass of 2.73 M Q . 
Compared with previous work [20], this value agrees 
within about 12%, and may be slightly overestimated. 
Other authors [H,[ll| often use a lower value of the coeffi- 
cient of nuclear incompressibility, namely K = 180 MeV, 
while we use K = 220 MeV. This significantly changes 
the maximum mass that the PNS can hold and thus, 
comparisons are more difficult. 

The model lsWW40 forms an AH at 1025.8 ms post- 
bounce, while the model 7rWW40 collapses earlier, and 
an AH is detected at 607.2 ms, for a maximum PNS mass 
of 2.40 Mq. We note that the addition of pions leads 
to larger changes in the properties at BH collapse for 
WWs40 than for u40 (see Tab.lHl). We may interpret this 
as follows: because the PNS contracts more when using 
LS220+7T EoS, the neutrinosphere is closer to the center 
and so, cooling is more effective. Because accretion rate is 
lower with the WWs40 progenitor, more time is needed to 
reach the PNS maximum mass. Hence, the more effective 
cooling lasts longer and differences between models using 
LS220+7T EoS and LS220 EoS are more pronounced. 

Fig. [9] has been translated with the respective bounce 
time of each model. Note that the bounce happens 
~ 12 ms earlier for model AWW40, as a small fraction 
of A begins to appear (Fa ~ 0.0005 at the center at 
bounce) before bounce. For this model BH collapse is 
triggered before the phase transition density is reached. 
Indeed, with progenitor u40 we reach a central density 
of 10 15 g.cm -3 at ~ 69 ms post-bounce, which triggers 
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FIG. 9: Central density as a function of time for model lsWW40 (dashed, green), 7rWW40 (plain, red) and AWW40 (dotted, 
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FIG. 10: Neutrino luminosity as a function of post-bounce 
time for model AWW40 (dashed, green), compared to ref- 
erence model lsWW40 (plain, red). Model AWW40 shows 
a very high peak luminosity and a more sustained neutrino 
burst. 



genitor that collapses to a BH. Note that, even without 
phase transition, we find a non-negligible fraction of A- 
hyperons (Ya ~ 0.15 at the center at the onset of BH col- 
lapse, at t = 271 ms post-bounce. See Fig. too, which 
shows Ya ~ 0.17 at the center in model Au40 before the 
phase transition). This is, as noted in Sec. IIV C 21 due 
to the fact that A-hyperons begin to appear before the 
phase transition at finite temperature, and their fraction 
increases with increasing temperature. 

Finally, Fig. [10] shows the total neutrino luminosity 
L„ (summed over all neutrino species) around the time 
of neutrino burst. Model AWW40 shows a high peak 
luminosity of 5.10 54 erg.s -1 , compared to the reference 
model lsWW40 that shows a maximum of 1.10 54 erg.s" 1 . 
In the model AWW40, the neutrino burst is also longer 
in time due to the rapid shrinking of the PNS, and the 
post-bounce luminosity stays higher until BH formation. 



V. CONCLUSION 



the onset of the phase transition. Here, as can be seen in 
Fig. [HI the collapse to the BH has already started when 
we reach a density of 10 15 g.cm~ 3 . This is due to the 
higher accretion rate of progenitor u40, which allows for 
higher central densities of the PNS. Hence, the LS220+A 
EoS does not induce a phase transition for every pro- 



We have presented simulations of stellar core-collapses 
to BH, comparing different finite-temperature EoSs with 
additional particles, namely pions and A-hyperons. As 
expected, additional degrees of freedom modify the EoS 
properties, in such a way that the collapse to a black hole 
occurs sooner after the bounce. 
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Our new EoSs are based on the Lattimer and Swesty 
EoS. Pions are added as a free gas, and A-hyperons are 
incorporated with the interactions adapted from Bal- 
berg and Gal [53[. This EoS is subject to a first-order 
phase transition driven by A-hyperons [36l . |54| which 
is described by a Gibbs construction. The LS220+7T 
EoS fulfills the gravitational mass constraint from the 
M = 1.97±0.04 M neutron star, and the LS220+A EoS 
is very slightly below with a maximum mass of 1.91 Mq. 
Compared to previous works, H.Shen et al. [34j j reported 
a maximum mass of 1.75 Mq for their EoS including 
A-hyperons, and Ishizuka et al. (30j reported a maxi- 
mum mass of 1.55, 1.63 and 1.65 M for their EoSs that 
include different parameterizations of hyperons, muons 
and pions. So, the EoSs presented in this paper are in 
better agreement with observational mass constraints. 

We have implemented a new leakage scheme. Com- 
pared with previous works, we keep track of the neutrino 
fractions more accurately with advection equations, and 
take into account neutrinos that are not trapped anymore 
because of the neutrinosphere moving inwards. Although 
our fluid rest frame source terms are approximated, we 
consistently transform them to the Eulerian frame, where 
the hydrodynamic equations are solved. With these re- 
finements, it is possible to follow a simulation during the 
collapse, bounce and post-bounce phases without chang- 
ing the approach. 

When using the LS220+A EoS, we have implemented 
for the first time the isoenergetic scattering off a A- 
hyperon. It is done in a simplified manner consistent 
with our other opacity sources, and we find that it shifts 
the occurrence of the phase transition by 6 ms and the 
PNS maximum baryonic mass is 0.04 Mq smaller with 
this opacity taken into account. This contribution is still 
very subdominant compared with scattering off nucleons. 

Our use of the leakage scheme during the collapse 
phase enabled us to see differences in the deleptonization, 
which resulted in a smaller density at bounce for model 
7ru40 and 7rWW40, compared to lsu40 and lsWW40, re- 
spectively. With both progenitors, we find, as expected, 
that the compression of the PNS is faster when using 
the LS220+7T EoS. The reduction of the PNS maximum 
baryonic mass is 0.06 M for model 7ru40 compared to 
model lsu40. A larger difference of 0.33 M Q is found when 
comparing models 7rWW40 and lsWW40. We interpret 
the larger difference as a result of the larger amount of 
time spent in the post-bounce cooling regime. Small dif- 
ferences between EoSs have thus more time to develop. 

With the u40 progenitor and the LS220+A EoS, a 
phase transition clearly appears. Contrary to the cold 
A EoS, at finite temperature A-hyperons begin to appear 
before the phase transition, and the fraction of thermally 
populated A-hyperons is significant (~ 0.17 at the center 
before the phase transition). The phase transition results 
in a sudden increase in density, comparable to the free 
fall during the initial collapse. 

In contrast to model Au40, with model AWW40 we do 
not observe the phase transition. We find that up to the 



onset of BH collapse, the density stays too low to trigger 
the phase transition. Consequently, the appearance of A- 
hyperons is entirely due to thermal effects. We conclude 
that our models with LS220+A EoS do not predict a 
phase transition for every progenitor. Indeed, only the 
progenitors with the highest accretion rates are able to 
reach the phase transition density before the onset of BH 
collapse. 

Apart from the lack of the phase transition, the 
WWs40 progenitor has the same qualitative behavior as 
u40. Our results with model lsWW40 are in agreement 
with previous work with the same progenitor and the 
same EoS [H. The LS220+?r EoS collapses to BH earlier 
than the LS220 EoS, admitting a lower PNS maximum 
mass. LS220+A EoS also leads to a more compressible 
PNS, and the effect is much more pronounced. The PNS 
maximum mass is found to be 2.00 Mq. 

Thus, we addressed the question of additional par- 
ticles in hot EoSs. It seems that adding pions or A- 
hyperons can significantly change the conditions of a 
core-collapse supernova, while still being in agreement 
with the M = 1.97 ± 0.04 M neutron star. Because we 
see differences already at bounce, and because the neu- 
trino luminosity is higher when adding A-hyperons, one 
could infer that additional particles may make a differ- 
ence in the explosion phase. This will have to be investi- 
gated with a more detailed neutrino treatment. 

The present EoS with A-hyperons is only marginally 
compatible with the mass of PSR J 1614-2230 and the 
question arises to which extent our results would be mod- 
ified by taking another EoS giving a higher maximum 
neutron star mass. As far as the same model is used, e.g. 
with one of the parameter sets from Ref. [35| , no qualita- 
tive changes are to be expected. Indeed, as can be seen 
from Ref. [35], the behavior of the hyperonic EoS with 
different parameterizations are very similar. More pro- 
nounced modifications are of course to be expected if a 
future observation gives an even higher neutron star mass 
and another model has to be used (remind that LS220 
EoS without additional particles has a maximum mass 
of only 2.06 Mq). The general effect of reduced time 
to black hole collapse in presence of additional particles 
seems, however, very robust, since we confirm the results 
from [30h32| | which have been obtained using an extended 
version of the HShen EoS, thus a different model for dense 
matter. 

Further studies on observational consequences of the 
appearance of hyperons or pions are yet to be done. In 
particular, the phase transition should produce copious 
amount of gravitational waves, as it has been shown by 
previous studies on "mini-collapses" [7(1 [77}. A future 
study shall use simulations of rotating stellar core col- 
lapse with pion and hyperon EoS, in order to infer the 
gravitational wave signal. On the other hand, our estima- 
tion of neutrino luminosity suggests that the PNS radial 
modes might be detectable, although it would be hard 
to unambiguously associate it to a phase transition to 
hyperonic matter, also because other phase transitions 
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(e.g. to quark matter [13]) are possible. So, it seems 
possible to detect evidences of the phase transition in an 
ideal case, combining observations of neutrinos and grav- 
itational waves from core-collapse to a BH, even without 
constraining the nature of the phase transition from these 
data. 

Our work on modern EoSs including pions and A- 
hyperons raises the question of the impact of the other 
additional particles that we did not consider. In the long 
run, it will be desirable to build an EoS with n, A but also 
other hyperons (£, S), and muons. This EoS would have 
to fulfill the mass constraint and consistently take into 
account interactions between particles. Coupling it to an 
accurate neutrino transport scheme, taking into account 
neutrino reactions with the new particles is the only way 
to have accurate quantitative results on the influence of 
additional particles. In the mean time, modern nuclear 
and astrophysical data should restrict more and more the 
set of compatible parameters, which will help building a 
realistic EoS. 



Newtonian limit that the Eulerian frame becomes the 
same as the coordinate frame. 

This justifies that the transformation from the LF to 
the EF is a Lorentz boost and contains no metric terms. 
Note, however, that our velocity Vi defined in Sec. Ill Al 
is written in the coordinate frame, and we need to trans- 
form it into the EF. 

To do so, we define the matrix of a general Lorentz 
transformation (tetrad transformation) v as 

g^P a ^Pp v = Va . (A2) 

Going from the coordinate frame velocity to the 
EF velocity is done by applying this transforma- 
tion matrix P to the 4-vector u M = (0,v l ). In 
the CFC, this simply results in multiplying Vi by 
diag(l/$ 2 , l/($ 2 r), l/($ 2 r sin6>)). 

Writing the boost for an arbitrary Eulerian velocity 
Vi,EF = {vr,ve,v<p) leads to 
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Appendix A: Explicit derivation of the 
hydrodynamic sources 

Neutrinos enter the fluid momentum and energy 
equations as a source term in the hydrodynamic 
Eqs. V^T^ = q a , see Eqs. flU EH- The quantities 
in Eq. ([M]) being defined in the Lagrangian frame (LF) , 
we further need to transform them in the Eulerian frame 
(EF), where the hydrodynamic equations, Eqs. (|8|). are 
solved. 

We adopt the same conventions as [40( and [58[ for the 
definitions of frames : The coordinate frame (that can 
be associated with the grid) is fixed, and its tetrad is 
used to define the metric as (see (icj |l 



9fiu = g(d^,d u 



(Al) 



Note also that the coordinate frame is not associated 
with a physical observer. The Eulerian frame is then 
defined as in Sec. Ill Al (the Eulerian observer moves or- 
thogonally to spacelike hypersurfaces) . It is an inertial 
frame, in which the basis vectors span the Minkowski 
metric rj a p = diag(— 1, 1, 1, 1), so that it is only in the 



<& F = K a p{v i > EF )<f LF (A3) 
and explicitly, the covariant components are 



q ,EF = W{Q E + V^ b Q(M) 



(A4) 



Ql.EF = v r WQ E 



q 2 = v B WQ E 



+ 



Qz,ef = V^WQe 



1 + {W-1) 



v r v 



V;V 



+ 



(W-l) 
(W-l) 



ViV 1 

v r v^ 

ViV 1 



!{M)l 
Q(M)2 

Q(M)3 (A5) 



(W-l) 



VgV 



Q(M)1 



l + (W-l) 



v e v 



v,v 



(W-l) 



VgV^ 



(W-l) 
(W-l) 



I'd, 



v 



ViV 
V^V 6 



Q(M)2 
Q(M)3 (A6) 

Q(M)1 



ViV 



Q(M)2 



1 + (W-1) 



<3(m) 3 (A7) 



It can easily be seen that Eqs. (|A4|). ([A5]) . (TAB]) . (|A7j) 
reduce to the usual spherically symmetric boost by im- 
posing vg = v^ — 0. Terms arising from the deviation 
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from spherical symmetry, while usually small, can be- 
come non-negligible, for instance in the case of a rapidly 
rotating core [58| , 

Finally, in Eq. ([8]), qa.EF enters as such in the source 
term of the conserved quantity E corresponding to the 
conservation of energy, and qj,EF enters as such in the 
source term of the conserved quantity Sj corresponding 
to the conservation of momentum. Note that because of 
the definition of Eq. ([8]). we recover there a multiplicative 
metric term, \J—g. 



b. Elastic scattering off a neutron 

Following Ruffert et al. , the opacity t s iyn) for the 
elastic scattering off a neutron is 



T s {vn) = C sn n n a a xYNNr=-, — r , 

m e c 2 F 3 (rj„) 



(B5) 



where n n is the neutron number density. The constant 
C sn takes the value 



Appendix B: Neutrino reactions formulae 

In this appendix, we do not use c = 1 for clarity. The 
temperature T is in units of energy. 



Csn - 24 



1 + 5|M 

gv 



(B6) 



1. Opacities 



c. Elastic scattering off a nucleus 



a. Elastic scattering off a proton 

The opacity T s (vp) is defined following Ruffert et 
al. [60| , with the cross section corresponding to the trans- 
port cross section from Burrows et al. [Hj]. 



T 

T s (yv) = C sp n p a 5- Y N N 



F 3 (Vu) 



(Bl) 



where n p is the proton number density, T is the tem- 
perature, r\ v is the degeneracy parameter, and uq is the 
weak interaction cross section, defined by 



4(m e c 2 G F ) s 
7r(7ic) 4 



(B2) 



with Gf being the Fermi constant, m e the electron 
mass, and F p the Fermi integral of order p defined by 
Eq. (1551 . The F 5 /F 3 term comes from phase space inte- 
gration and corresponds to the energy opacities in [60| . 
which are the ones we find in better agreement with more 
detailed calculations. We also include the following ap- 
proximation for the nucleon-nucleon degeneracy factor 
Ynn, 



We follow Rampp and Janka 78] here, and define the 
opacity t s {uA) as 



t s {vA) = — -n A - 



16m e c 2 F 3 (r) v ) 
x ( (C A - C v ) + (2 - C A - C v ) 



2Z - A 



A 



x tt- 1+ (Y*> e ~ a 'V( e ))< 



(B7) 



where n A is the heavy nuclei number density, and 
C A = \. The ion screening term (S(e))i on is taken 
into account as in Horowitz [79( and is defined as in 
Bruenn and Mezzacappa (80l | and Rampp and Janka [78[ , 
cxplicitely 



Vb = 



5{hc} 



r (1.07A) 2 / 3 < 



(B8) 



with A the mass number of the average heavy nucleus, 
and (e u ) the mean neutrino energy, defined by Eq. 1371 



Y, 



Y, 



NN 



1 



max(^,0) 



(B3) 



where Y/y is the nucleon fraction, [a^ is the nucleon 
chemical potential and T is the temperature. Finally, 
the constant term C sp takes the value 



C 



sp 



1 

24 



4(C V 



1) 



9A 

gv 



(B4) 



where gv and g A are the vector and axial coupling 
constants of the weak interaction, and Cy = 5 +2 sin 2 9\y 
with 9w the Weinberg angle. 



d. Absorption of av^ by a neutron 
The opacity T a (v e n), from [601 ] . for this reaction is 



T a (y e n) = — 



1 + 3| ^ 

gv 



F 5 (y,J 



1 + exp 77, 



colKd] 
F 5 (Vu, 



T 



-1 



Fii-q, 



(B9) 



with rj e the degeneracy parameter of the electrons and 
\Vud\ the CKM matrix element. 



19 



e. Absorption of a v e by a proton 
The opacity r a (P e p), from (60l |. for this reaction is 



Ta(y & V) 



1 + 3 



0o|Kd| 



T 



^5 fee) 
^3 fee) 



1 + CXp -77 e — 



f 5 fee) 
Fi fee) 



n -1 



(BIO) 



These integrated rates are computationally expensive, 
and only require the knowledge of the EoS to be com- 
puted. Therefore, they can be tabulated as a function 
of p,T,Y e (at each point of the EoS), and the effective 
neutrino chemical potential fJ>u,eff- Access to the tables 
is done by a quadrilinear interpolation. 

j is the creation term, also taken from Bruenn [64l |. 



b. Electron capture on free protons 



f. Elastic scattering off a A-hyperon 

Following the same approximations as for the transport 
cross section for coherent scattering off neutrons in [12|, 
IH H we obtain, 



CTA 



cr /l + 5cj 



A A 



6 



(Bll) 



The constant ca a is the A axial coupling constant to the 
neutral current and is assumed to take its value given by 
flavor SU(3) symmetry, c Aa = -0.73 [HIH 2 . 
Explicitely, the opacity r s (vA) for this reaction is 



t s (vA) = n A <j A 



T 



(B12) 



where n A is the A-hyperon number density. Note that 
we do not take into account a degeneracy factor in this 
case. 



2. Neutrino creation terms 

a. Electron and positron captures 

To take into account electron and positron captures, 
we integrate the following expression (see (63 |) 



47T 



n b (hcf 



/>oo 

/ E 2 j(l-f)dE, 
Jo 



(B13) 



where £ = E ec ,Sp C ; / is the neutrino distribution 
function, which we assume is a Fermi-Dirac, because no 
deviation from equilibrium is computed in the leakage 
scheme. Note that we do not take into account the ab- 
sorption term, and rather set /3-equilibrium at a given 
density. 



The total electron capture rate is the sum of the elec- 
tron capture rate on free protons and the electron capture 
rate on nuclei. For electron capture on free protons, j eCtP 
takes the value 



^jlKdlVfffv + 3<?i)/ e (£ + Q) 



x {E + Qf 



1 - 



m e c 
E + Q 



1/2 



(B14) 



where Q is the mass difference between neutrons and 
protons, and r\ pn is such that 



, d 3 p 
'(2tt) 3 



f n (E)(l-f p (E)) 



(B15) 



with / e , f n and f p the distribution functions of the 
electrons, neutrons and protons, respectively, all taken 
to be Fermi-Dirac distribution functions. 



c. Electron capture on nuclei 
For electron capture on nuclei, j ec . n takes the value 



a] 



Je 



Tr(hc) 4 
x (E + Q') 2 



n A \V ud \ 2 g 2 A N p (Z)N h (N)f e (E + Q') 
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with the approximation Q' = ji n — fj, p + A, /j, n and /i p 
being the neutron and proton chemical potential, respec- 
tively, and taking the constant value A = 3 MeV. 

Np(Z) and Nh(N) are also taken following Bruenn f&il ]. 



0, Z < 20 

N p {Z) = <| Z- 20, 20 < Z < 28 
8, Z > 28 



(B17) 



2 Note that the strangeness content As of the nucleon has been 
neglected when deriving this value, see [8l| . The value of As ~ 
—0.1 is subject to very large uncertainties and taking it into 
account would only very slightly modify our results. 



6, N < 34 

N h (N) = { 40 - N, 34 < TV < 40 
0, N > 40 



(B18) 
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We start the integration Eq. (|B13|) at the threshold 
E = m e c 2 — Q' if this value is positive. 



d. Positron capture on free neutrons 

For positron capture on free neutrons, j pc takes the 
value 



The constant C ee takes the value C ee = ((CV — Ca) 2 + 
(CV + C A ) 2 )/36 for v e and D e , and C ee = ((CV - C A f + 
(C v +C A - 2) 2 )/9 for v x . 



f. Plasmon decay 



G 



Tr(hc) 4 
x {E-Qf 



1 l Kd|V P (<7y + 3.gi)/ e+ (S-0) 



m e c 
E-Q 



-, 1/2 



(B19) 



where f e + is again a Fermi-Dirac distribution function. 
rj np can be inferred from the definition of rj pn by inter- 
changing n and p. We start the integration Eq. (|B13p at 
the threshold E = Q + m e c 2 . 



e. Electron positron pair annihilation 
Following Ruffert et at. [U, 



J ee — ^ee 



CTQC 



(m e c 2 ) 2 e 
with the definition 



-e e +(l-A>ee(l-/ Ps >ee , (B20) 



Following Ruffert et al. [6C 



S p i — C 2 



a c 



ji8 



3a R 



m e c 2 ) 2 (he) 6 



7 6 e-^(l+ 7 )(l-/, i ) 7 (l-/ p ,) 7 , (B23) 



by 



where a em is the fine structure constant, 7 is defined 



7-2 



3tt 



r (7r 2 + 3r? e 2 



(B24) 



and the average Fermi distributions are 



87T 4 / 

and the average Fermi distributions 



(B21) 



v 1 - fui)*, = 1 +cxp 



- 1 



7 



1 + exp 



F 4 (77 e ) F 4 (-77e) 



2F 3 (r? e ) 2F 3 (-r; e ) 



2(1+7) 

(B25) 

Note here that i/ x are taken into account with a de- 
generacy of 4 and a vanishing chemical potential. The 
constant C p i takes the value C„i = Cy for v e and v ei and 



(B22) C p i = C V -1 for ^ 
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